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Abstract 

For a system at given temperature, with energy known as a function of a set of variables, we 
obtain the thermal fluctuation of the evolution of the variables by replacing the phase-space with 
a lattice and invoking the principle of detailed balance. Besides its simplicity, the asset of this 
method is that it enables us to obtain the Langevin equation when the phase-space is anisotropic 
and when the system is described by means of curvilinear coordinates. As an illustration, we apply 
our results to the Kramer- Watts- Tobin equation in superconductivity. The choice between the Ito 
and the Stratonovich procedures is discussed. 

PACS numbers: 05.10.Gg, 05.40.-a, 74.40. +k 



1 



I. INTRODUCTION 



In this article we deal with what we call "purely diffusive systems" in contact with a 
thermal bath. By this we mean systems with a state described by microscopic variables 
Xi, . . .xtv, with energy E{xi, . . .xn), which in the absence of thermal fluctuations are ex- 
pected to follow an evolution equation 

dxj dE 

— = -r,(x„...x.)— , (1) 

where t is the time and the positive coefficients Tj are determined by the dynamics of the 
system and its interaction with its environment. dE/dxj may be regarded as a driving force 
and Tj as a compliance coefficient. The Langevin approach tells us that the influence of 
thermal fluctuations can be taken into account by adding a fluctuating quantity at the right 
hand side of the evolution equation; this fluctuating quantity is called the "Langevin term." 

The paradigm of a purely diffusive system is a particle that undergoes Brownian motion. 
In this case the variable is its momentum, E is its kinetic energy and F is the Stokes 
coefficient. Paul Langevin dealt with this problemi and determined the variance of the 
Langevin term by invoking the theorem of the equipartition of the kinetic energy among 
the various degrees of freedom of a system in thermal equilibrium. Gillespie^ notes that it 
is not obvious that the influence of fluctuations can be separated as an additive term with 
zero average; regarding the velocity evolution as a Markovian process and assuming that 
the "stepping functions" (will be deflned in the following section) are linear functions of 
velocity, it is shown that this separation indeed occurs. Katayama and Terauti^ used the 
Langevin equation to study Brownian motion of a single particle under steady plane shear 
flow. Balescu^ introduces a Langevin equation in a model for the description of a plasma. 
Bringuier- discusses the difficulties encountered when applying the Langevin approach to 
the Hall effect. The relation between the Langevin and the Klein-Kramers approaches is 
discussed in Ref. [6. 

During the century that has elapsed since Langevin's paper, his approach has been ex- 
tended to wide classes of problems in Physics^"^'' and the Langevin term is determined by 
means of the fluctuation-dissipation theorem.—""— In this paper we will evaluate the distri- 
bution of the Langevin term by means of the principle of detailed balance.— 

The tools for handling the problems considered in this paper may be found in the literature 



on stochastic differential equations (e.g. Refs. 



and many of the questions raised here 



may be avoided by switching to the Fokker-Planck equation; this article is addressed to 
those readers, pressumably physicists, who prefer a more intuitive approach. 



II. OUR METHOD 
A. ID case 

We consider a one- dimensional system with microscopic state determined by the variable 
X. We discretize x and assume that its possible values are Xi = ii, where £ is a "lattice 
constant." We denote by = E{xi+i) — E{xi) the energy difference between consecutive 
lattice points. We assume now that for a short period of time 6t the probability of passing 
from Xi to Xi±i is given by Wi±6t, where Wi+ and Wi- are "stepping functions" that stand 
for the transition probability rates. 

The principle of detailed balance asserts that in thermal equilibrium the probability for a 
transition from i to i + 1 equals that for a transition in the opposite direction, i.e., denoting 
by the equilibrium probability for the value x = Xi, Pl'^Wi^^ = P^^-^^Wi+i^^. Since 
peq^peq^ _ ^ei/ksT ^ where ks is the Boltzmann constant and T is the temperature. 



It should be noted that not every system obeys detailed balance. Denoting by Pj the 
probability for x = xi (not necessarily for equilibrium), PiWi^+ — Pi+iWi+i- stands for the 
probability current. Detailed balance requires that this current vanishes, whereas in order 
to mantain a stationary state it is sufficient that the divergence of the current vanishes. If 
probability currents are present in equilibrium, it follows that there are driving forces which 
cannot be expressed as the gradient of the energy as in Eq. ([T]) (such as the magnetic force 
on a charged particle). Therefore, systems that do not obey detailed balance are beyond the 
scope of this article. 

In order to obtain more symmetric expressions, we write Wi^+ = WiXi, Wi- = Wi/Xi 
and, taking i sufficiently small so that quantities of order 0(£^) can be dropped, we write 
Wi±i = Wi±Wi and Ai±i = Aj ± where w'Jwi and A^ are at most of order 0{i). With this 
notation and approximation, Eq. ([2]) becomes 



(2) 



Wi + w[ 



(3) 
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Similarly, requiring detailed balance between the sites i and i — 1 we obtain 

{w. - w'^X, - K)e-/'-'^ = ^ , (4) 

where we have exchanged sides in the equation and neglected the 0(£^) difference ej_i — e^. 
Multiplying the equations ([3]) and (jl]) and neglecting the 0(£^) term {w'JwiY we obtain 

hence is of order 0(£^) and will be dropped. Aj can now be obtained from Eq. Q; keeping 
terms of order 0{i) and making use of the definitions of w'^ and e,, it becomes 

Xi = ^J{1 + w'Jwi)e--^I^BT = 1 + w'J{2w,) - ti/{2kBT) 

= 1 + {i/2)d[\og{w)]/dx - [e/{2kBT)]dE/dx , (6) 

where we have defined a smooth function w such that w{xi) = Wi. 

Let us denote by 6x the increment of the variable x during the period of time 6t. For 
sufficiently small 6t we can neglect multiple transitions and the possible values of 6x are 
and ±i. The average value of 6x will be 

(Sx) = iwiSt{Xi - 1/A,) = e'^w6t[d{\ogw)/dx - {l/{kBT))dE/dx] , (7) 

where in the last step we have neglected higher orders of i and have dropped the index i. 
Similarly, the variance of 6x will be 

{{SxY) = fwi6t{Xi + 1/Xi) = 2fw6t , (8) 

where besides dropping terms that are of higher order in i we have used the fact that, for 
small St, (Sx)'^ < {{S^Y)- 

We now get rid of the unphysical lattice by defining T{x) = f"w{x) / {kBT). With this 
notation Eqs. ([7]) and ([8]) become 

{5x) = T{x)6t[kBTd(logT)/dx - dE/dx] (9) 

and 

{{6xf) = 2kBTT{x)5t . (10) 
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Finally, we consider a lapse of time r which is very short compared with the relaxation time, 
but very long compared with 5t. By the central limit theorem,— li^ii^ the increment of x 
[which is the sum of many increments described by Eq. ^] will be 

Aa; = T{x) [ksTdilog T) jdx - dE/dx]T + rj , (11) 

where 1] is a fluctuating term with average 0, variance 2kBTT{x)T and Gaussian distribution. 
7] is the Langevin distribution we were looking for. There is still a subtle question concerning 
the precise value of x at which T{x) has to be evaluated; this issue is considered in Appendix 

m 

Let us now compare the nonfluctuating part of Eq. (ITT!) with Eq. ([1]). If F is independent 
of X, Eq. ([1]) is recovered; otherwise, there is also a drift term kBTd{\ogV) / dx that pushes x 
towards values where F is larger. The drift term can be absorbed into Eq. ([T]) if we replace 
the energy E hj G = E — /c^TlogF. The term /c^TlogF may be regarded as a sort of 
chemical potential, where F plays the role of the activity. 

B. Multivariable system 

We consider now a system with variables xi, . . .xn- For each of the variables we can 
repeat the analysis of the previous section and Eq. ffTTj) generahzes to 

Axj = ~Tj{dGj/dxj)T + rij , (12) 

with Gj = E — /c^TlogFj and rjj is a Langevin function with average 0, variance 2kBTTjT 
and Gaussian distribution. 

We might also be interested in the evolution of other variables rather than those in the 
set Xi, . . . xjv- This problem is considered in Appendix [Bl 

C. Curvilinear coordinates 

Now the volume in phase space is proportional to the Jacobian of the coordinates; there- 
fore, different lattice points in the discretized phase space may represent different volumes. 
As a consequence, P'^'^ has to be multiplied by this Jacobian and equations like Eq. ([2]) have 
to be modified accordingly. Before we deal with the general case, let us consider the case of 
polar coordinates. 
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1. Polar coordinates 



Let the coordinates be r, with volume element rdrd^p. For the variable the analysis 
remains unchanged, but for r we have P^'^/Pl^^ = e^^^'^^'^r / (r + i) and equations like Eq. 
have to be replaced with equations like 

= Wi^+e'^/'^^r/{r + i) . (13) 

Following the steps of Sec. Ill At instead of Eq. fllip we now obtain 

Ar = TrdlksT log(rr^) - E] jdr t + rjr , (14) 

where 1]^ has average 0, variance 2kBTTrT and Gaussian distribution. The replacement of 
the term logF^ with logrF^ implies a drift towards larger values of r. 

It is temptingi^ to attribute the drift towards larger values of r to the fluctuations of (p: 
if the system moves in phase space by the amount r5ip perpendicular to the radial direction, 
the new value of r would be ^Jr"^ + {r5i^Y ~ ^"[1 + ('^V^)^/2]- After a lapse of time r this effect 
would contribute an increment of r by the amount r{{Aip)'^ /2) = rksTT^T. Comparison 
of this result with Eq. ( IT^ indicates that this interpretation would be consistent with the 
principle of detailed balance only if F^ = r^F^. 



2. General case 



Let the coordinates be f i, . . . vn, with volume element J(f i, . . . VN)dvi, . . . dv^- Then, 
when dealing with the transitions of the variable Vj, Eq. (|T3|l generalizes to 

Wi+i^_ = Wi^+e'^/''^^J/{J + idJ/dvj) . (15) 



Following the steps of Sec. Ill At Aj has an additional term [i/(2J)]dJ/dvj and Eq. ([T] 
generalizes to 

Avj = Tjd[kBT \og{JTj) - E]/dvj t + r]j , (16) 
where rjj has average 0, variance 2kBTTjT and Gaussian distribution. An analogous result 



for macroscopic variables was obtained in Ref. 



20 



A Mathematica-program that illustrates the use of this result is provided in Appendix O 



6 



III. APPLICATION A MODEL FOR SUPERCONDUCTIVITY 



One of the most useful models in the study of dynamic properties of superconductors is 
the time-dependent Ginzburg-Landau model.— In this model the microstate of a super- 
conductor is described by a complex field ipi^x^y^z), such that ?/, 2;) p is proportional 
to the density of superconducting electrons at position {x,y,z). Knowledge of the field ip 
enables us to evaluate several measurable quantities, such as the supercurrent density. Since 
in this model the variable of the problem is itself a field, it may provide an example in which 
the Langevin approach appears to be more practical than the Fokker-Planck equation, since 
the latter is a partial differential equation in a space with infinitely many dimensions. 

In most cases, the time- dependent Ginzburg-Landau model is justified only for temper- 
atures very close to the transition temperature. The model was generalized by Kramer 
and Watts- Tobin;^'^ this generalized model is expected to be valid as long as there is local 
equilibrium. In order to focus on the aspects that we want to illustrate, we deal here with 
a simplified situation of the Kramer-Watts- Tobin model. We consider a uniform ID su- 
perconductor with periodic boundary conditions and ignore the electromagnetic field. We 
discretize the system by dividing it into N segments of equal length and denote by ipj the 
value of ijj at segment j. With appropriate normalizations, the energy of the system is given 
by 

^ 1 

where ^ is a constant that depends on the material, the temperature, and the length of each 
segment. If fluctuations are ignored, the evolution of ipj is given by^ 

fd^, + t^^^ = (1 - 1^/)^. + + V^.-i - 2^.) ' (18) 



where u and 7 are additional positive constants of the model. 

Let us first consider the case 7^(|V'jP) ^ 1, so that the left hand side in Eq. ( ITSj) can be 
approximated by udipj/dt. In this case it is convenient to express ipj in Cartesian form, ipj = 
Xj+iyj and the energy becomes X]iLi[~(^i+l/i) + (^j+?/j)^/2+^^(a;j-Xj_i)^+^^(?/j-yj_i)^]. 
Performing the derivatives and separating real and imaginary parts, Eq. (|T8|l takes the form 

^ = -r.|^ , (19) 

dt dxj ' ^ ' 
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with = l/{2u), and an analogous equation is obtained for yj. Since u is constant, we can 
apply the result f lTTj) with no drift term and conclude that after time r the fluctuating part 
of the increment of Xj will have a variance ksTT/u. 

Let us now consider the general situation. In this case it is convenient to express tpj in 
polar form, ipj = rje^'^\ the energy becomes X]^i[~''"| + ^j/2 + ^^('^| + ^|-i~2rjrj_i cos((pj — 
ipj^i))] and the expression in brackets at the left hand side of Eq. (fT8|) becomes [{drj / dt){l + 
7^r|) + irjd(pj/dt]e^'^K Multiplying Eq. f|T8|) by e"^'^^ and taking the real part we obtain the 
evolution of r^, 

dvj 1 dE , , 

-IT = / TT- 5 20 



taking the imaginary part gives the evolution of ipj 



dt c?99,- 



(21) 



These equations are in the form of Eq. ([T]), with = l/(2n^/l + 7^r|) and F 



1 + 'j'^r'j / {2urJ) . In the extreme case 'yrj ^ 1, F^ ~ l/(2u'~frj) and F^ ^ '~f/(2urj). 
Since in this limiting situation rjF^ does not depend on rj and r^F^ does not depend on <fj, 
the correction terms log(rF) are not required and the formalism developed in Sec. Ill CTl can 
be applied with no drift terms. For general 7, thermal fluctuations add a drift to Eq. f l2U]) 
and lead to 



1 + 'y'^r'^ J dr 



T + r]r 



ksT dE 



T + Vr. (22) 



where 77^ is the usual Langevin term with variance 2kBTTrT. 

It should be emphasized that Eq. f l22|) (including the drift) is an extension of Eq. f l20|) 
and not a modification of it. As an illustration of this statement, let us focus on the case 
7 = already considered in Eq. ( fTOl) . In this case the drift term becomes ksT/vj 7^ and 
Tr becomes 1/(2m), i.e. F^ = T^- Moreover, let us for a moment leave the KWT model aside 
and consider the toy model E = J2f=i — i-^-' J^^t have udipj/dt = —ipj 

instead of Eq. f lTSj) and udrj/dt = —Vj instead of Eq. fl20|) . It follows that if the drift term 
ksT/rj were not present, the evolution equations for Xj, yj and rj (including fluctuations) 
would all become identical and we would therefore have (r|) = = whereas the 
true relationship is (r|) = (x| + yp = 2{xp. 
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IV. CONCLUSION 



We have developed a simple method that enables us to derive the Langevin and drift terms 
for systems with random-walk type evolution. The cocepts of probability theory that we 
have invoked are intuitive, elementary and in the "language" used by undergraduate Physics 
textbooks. This method is particularly useful when the dynamics leads to an anisotropic 
phase space or when the evolution is naturally expessed in curvilinear coordinates. 
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Appendix A: Ito or Stratonovich? 

We first note that both terms in Eq. (11 II) are not of comparable sizes. The first is of order 
0(r), whereas rj is mean that rj is more important, since it tends to 

cancel in the long run, whereas the first term persists. 

The next question concerns the precise value of x at which Eq. ( ITT]) should be evaluated. 
Should it be the initial value, x = Xjn, the final value x = Xi^ + Ax, or some intermediate 
value? For the first term this question is irrelevant, since the change in the value of this 
term is 0{tAx) and its contribution vanishes in the limit r — t- (even after noting that 
decreasing the time interval by some factor increases the number of intervals by the same 
factor). This is not the case for the choice of the x-dependent variance of r]. 

The Ito procedure evaluates the variance at x = x-m- f] has a symmetric distribution and 
(77) = 0. The Stratonovich procedure evaluates the variance at the middle of the interval, 
X = Xin+Ax/2. In order to distinguish between the two procedures, we denote the respective 
random terms by 771 and In order to fix ideas, let us explore the case that T{x) is an 
increasing function of x. Since the variance is evaluated at x ~ Xin + ris/2, it means that 
if r^g > (respectively r^g < 0), then rjs/rji > 1 (respectively rjs/rji < 1). Qualitatively, 
this means that the distribution of rj^ will have a longer tail than that of rji in the positive 
direction and a shorter tail in the negative direction. 
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In more quantitative terms, we can establish a one to one correspondence between the 
values of r/s and those of r/i, such that 

Vs ~ [r(xin + r/s/2)/r(a;in)]'/V ^ {l + T'r]s/2TY/\ ^ {l+T'r]i/AT)T]i = rn + T'rjf/AT , (Al) 

where we have written F as shorthand for r(xin) and F' for dT/dx at x = Xin- The approx- 
imations in sequel ( lAip are Ax ~ r^s, expansions to first order in Ax, and risf]! ~ rji. The 
important consequence is that for sufficiently small r we have 

(r/s) = (r/i) + (r74r) {r^f) = kBTV'r/2 . (A2) 

This term is not negligible and does not cancel in the long run. 

In a didactic article, van Kampen^ explains how to translate between the Ito and the 
Stratonovich procedures. He advocates the use of a master equation rather than a Langevin 
approach, so that the Ito-Stratonovich dilemma never arises. Langon et al. performed an 
experiment in which colloidal particles diffuse in a medium with position-dependend diffusion 
coefficient, so that (r/s) 7^ {rji). In their case they found that r(x) has to be evaluated at 

In order to judge what is the appropriate procedure for Eq. f|TT]) in our case, we will 
evaluate (77). Since the difference between both procedures depends on w{x) and not on A, 
we are free to take A = 1 (we may imagine that E — /c^TlogF is constant) and are left 
with a random walk problem in which each step is equally probable for both directions. Let 
the system be initially at x = Xm; by definition, the probability distribution for 77 is the 
probability distribution for x — Xin after time r. 

Let Pi be the probability to find the system at x = if at some moment. The change 
of probability after time 5t will be 5Pi = {wi-iPi-i — 2wiPi + Wi^iPi+i)6t. Or, defining 

Pi = WiPi, 

5pi = Wi{pi_i - 2pi + pi+i)5t . (A3) 

In order to use the central limit theorem, we require 5t/r — )■ and we pass to a continuous 
model. Equation (]A3|) becomes the diffusion equation 

In order to pass from Eq. (1A3|) to Eq. (lA4p we have expanded p(x) to order 0(£^) and used 
the definition of F. Equation (1A4p may be regarded as our master equation. 
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FIG. 1: Deviation of the probability density from the Ito distribution. The unit of length can be 
chosen arbitrarily. The unit of time equals the square of the length unit divided by kBTT{xi^). 
The probability density is in units of d\ogT/dx. 

Let us take the initial value 0) = w{x\-a)5{x—x\-a). If in Eq. (1A4I) r(a;) were substituted 
by the constant r(xin), the solution of the diffusion equation would be 



(0)/ ,x ""^(^^in) 

'{x,t) = — exp 



i^X XiYl) 



2 1 



AkBTT{x,^)t^ 



(A5) 



2^'KkBTr{x,^)t 

which is just the Ito distribution multiplied by w{xia). We now deal with Eq. (IA4I) by means 
of the approximation 

do d'^D d'^o 92 ^(0) 

^ = kBT[V{x,^) + r(x) - r(xi,)]^ ^ kBTV{x,^)^ + kBTV ■ (x - a^in)^ . (A6) 

This nonhomogeneous equation is solved using the Green function of the diffusion 
equation.— We obtain that the deviation of the probability distribution from the Ito dis- 
tribution is AV{x,t) = ((ilogr/(ix)(— 60 + 0^) exp(— 0^/4)/(16yvr), where = (x — 
Xin) / >ykBTT{xin)T and the derivative is evaluated at x = Xin. A plot of this deviation 
is shown in Fig. [1] [for any value of r for which the approximations in Eq. (]A6|) are justified]. 

The deviation AV enhances the positive tail of the distribution of rji and hinders the 
negative tail, as in the case of 1]^. However, the key feature is that {r]) is not affected by AV; 
therefore, rji has the appropriate distribution for our problem and the Ito procedure should 
be used. Also {rj'^) is not affected by AV. It is interesting to note that, no matter how small 
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r is, there are always values of x where AV remains finite; however, the statistical weight 
of AV does become negligible in the limit of small r. 



Appendix B: Linear transformation of variables 

Let xi,...xn be the original set of variables for which the evolution of the system is 
known. Let us restrict ourselves to cases in which the coefficients are independent of the 
coordinates, so that Gj = E. Let Ui, . . .un he a new set of variables, defined by means of a 
linear transformation 

N 

where Cij are the elements of a constant orthogonal matrix C; the purpose of the factor 
yT^ is to compensate for the anisotropy of phase space. 
Substituting Eq. ([12]) into Eq. ([BT]) we obtain 

N op 

A«, = -5^c,,Vr;^^ + 6, (B2) 

where = Ylf=i ^ijVj/ is a fiuctuating term with zero average and Gaussian distribu- 
tion. Its variance is 

TV N 

J2 CUv])/^j = ^ksTr J2 Cl = 2kBTT , (B3) 
i=i 3=1 

where we have used the property that rjj and rjji are not correlated for j ^ j' and orthogo- 
nality of C. For i ^ i' we have 

N N 

where we have used again {rjjTjji) = 2kBTTT j5jji and orthogonality of C. 

Since Xj = a/F" Y.iLi ^jl'^'i = a/Tj J2iLi CijUi, dxj/dui = ^/r~Cij. Therefore, 

OUi ^ OUiOXj ^ OXj 



Comparing Eq. (1B5I) and Eq. ( 1B2I) we finally obtain 



Am, = -— r + ^i. (B6) 

OUi 
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Appendix C: Numeric Test 



In this Appendix we evaluate the statistical average of several quantities for a 2D system, 
using polar coordinates. In the following example we take an energy in which the Cartesian 
coordinates separate, namely, E = x"^ + 2y'^; the statistical average of any function of x will 



be (fix)) = J^^f{x)exp[-x'^/{kBT)]dx/ J^^exp[-x'^ /{kBT)]dx and, similarly, {g{y)) = 
JZ,g{y)exp[-2yy{kBT)]dy/JZ,exp[~2yy{kBT)]dy and {f{x)g{y)) = {f{x)){g{y)). The 



following is a Mathematica-program that evaluates several averages of this kind. 
Clear [r, phi] ; (* These are the polar coordinates r,(p *) 

energy = r"2 (1 + Sin [phi] "2); (* This is the energy in polar coordinates *) 
kT = 2; (* This is ksT *) 

Gammartau = 5.*"-4 r; Gammaphitau = 5.*~-5/r"2; (* These are T^t and T^t; 
only the products of these quantities appear in each step; statistical averages should be 
independent of these dynamical functions, which have been chosen arbitrarily *) 
stepr = Simplify [Gammartau D[kT Log[r Gammartau] - energy, r]] ; (* This is 
Ar as in Eq. f lT^ without r]r *) 

stepphi = Simplify [Gammaphitau D[kT Log[r Gammaphitau] - energy, phi]]; (* 
This is A(p without t]^ *) 

stdr = Simplify [Sqrt [2 kT Gammartau], Assumptions -> r > 0] ; stdphi = 

Simplify [Sqrt [2 kT Gammaphitau], Assumptions -> r > 0] ; (* These are 
the standard deviations of rjr and t]^ *) 

TABLE I: Statistical averages of several quantities, for two temperatures, "theory" refers to 
ensemble averages and "numeric" to average over steps, as obtained from the evolution predicted 



by Sec. HTCTl 



{xy) {\x\) (x2) {\y\) (y2) 



theory, ksT = \ 



0.00 0.564 0.500 0.399 0.250 0.225 



numeric, ksT = 1 



-0.01 0.561 0.493 0.402 0.253 0.227 



theory, ksT = 2 



0.00 0.798 1.000 0.564 0.500 0.450 



numeric, ksT = 2 



-0.02 0.813 1.028 0.557 0.490 0.456 
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r = 0.1; phi = 0.2; (* This is the initial microstate; the statistical averages should 
be independent of it *) 

Nrelax = 5 10"6; (* Number of steps during which the system "forgets" the initial 
state and relaxes from it to a "typical" microstate *) 

Do [r = r + stepr + RandomReal [NormalDistribution[0, stdr]]; (* r evolves 
during a step, according to Eq. (IT^ . stdr is evaluated at the beginning of the step, 
according to the Ito procedure *) 

If [r < 0, r = -r; phi = phi - Pi] ; (* If r becomes negative, r and <f are 
redefined *) 

phi = phi + stepphi + RandomReal [NormalDistribution[0, stdphi]], (* (p 
evolves during a step, according to Sec. Ill C 1[ does not depend on (p, so that the stage 
at which stdphi is evaluated is not crucial *) 

{i, Nrelax}] ; (* Length of the loop *) 
Naverage = 15 10"6; (* Number of steps during which averages will be evaluated *) 
sx = 0; sxx = 0; sy = 0; syy = 0; sxy = 0; sabs = 0; (* Initialization of the 
variables that will be used for the evaluation of cumulative sums, from which the averages 
will be obtained; the following loop is identical to the one above, except that now we keep 
track of these sums *) 

Do [r = r + stepr + RandomReal [NormalDistribution[0, stdr]]; 
If [r < 0, r = -r; phi = phi - Pi] ; 

phi = phi + stepphi + RandomReal [NormalDistribution[0, stdphi]]; 
X = r Cos [phi] ; y = r Sin [phi] ; sxy = sxy + x y; x = Abs [x] ; 
y = Abs [y] ; sx = sx + x; sxx = sxx + x"2; sy = sy + y; 
syy ~ syy y"2; sabs = sabs + X y, {i, Naverage}] ; 
Print ["<xy>=" , sxy/Naverage] ; Print ["< I x I >=" , sx/Naverage] ; 
Print ["<x"2>=" , sxx/Naverage] ; Print ["< I y I >=" , sy/Naverage] ; 
Print ["<y"2>=" , syy/Naverage] ; Print ["< I xy I >=" , sabs/Naverage] ; (* The aver- 
ages that we decided to evaluate are printed *) 
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In Table H] we compare the results obtained by this program with the expected statistical 
averages. 
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